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ABSTRACT 


A parametric studv was conducted to assess the effect of thermal residual stresses 
on the stress-strain response of a fiber-reinforced metal-matrix composite. The material 
chosen for investigation was the SiC-whisker reinforced Al 6061 system. The eflects of 
fiber volume fraction, fiber aspect ratio and fiber spacing were analyzed within the 
framework of axisymmetric finite-element models to determine the overall constitutive 
response of the composite as well as to solve for local field quantities in the fiber and 
matrix. The composite was modeled as a periodic array of cylindrical fibers, laterally 
aligned in one model and staggered in the other. Perfect interfacial bonding and com- 
plete fiber alignment with the tensile axis were assumed. The results indicated that (1) 
composite stiffness, yield strength and work-hardening rate increased with increasing 
volume fraction and fiber aspect ratio and (2) variations in fiber spacing primarily affect 
work-hardening rate and have negligible effect on composite stiffness. lt was found that 
the presence of residual stresses affected the stress-strain behavior of the composite by 
imfluencing the load transfer characteristics between the matrix and fiber as well as the 


initiation and growth of plastic deformation in the matrix. 
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I. INTRODUCTION 


Metal matrix composites are materials that combine a relatively soft metallic matrix 
with a high-strength, high-stiffness reinforcement material that is tvpicallv a ceramic, 
intermetallic, metalloid or metal. These materials are characterized either bv the nature 
or the morphology of the reinforcement. Continuous fiber composites consist of fibers, 
often in a laminar construction, spanning the entire length of the matrix material. Dis- 
continuous fiber composites consist of short fibers or whiskers arranged in a random 
fashion in the matrix. The degree of randomness as well as the length of the fibers can 
be controlled, to an extent, through processing. Other metal matrix composite systen1s 
such as particulate, flake and foam composites derive their names [rom the morphology 
of the reinforcing material. 

Wrought aluminum alloys, such as 2024, 6061 and 7075, reinforced with boron, 
graphite, alumina or silicon carbide have been used extensively in applications as engine 
Or structural components. Some of the improvements in mechanical properties over the 
unreinforced matrix include increased elastic modulus, tensile strength, wear resistance, 
fatigue resistance, strength-to-weight ratio and thermal conductivity. 

In order for optimal utilization of metal matrix composites in industry, detailed 
Knowledge of the factors affecting their thermal, electrical and mechanical properties 


must be known. 


A. THERMAL RESIDUAL STRESSES IN METAL-MATRIX COMPOSITES 
Thermal residual stresses are generated in metal matrix composites during fabri- 

cation due to the coefficient of thermal expansion (CTE) mismatch between the matrix 

and reinforcing material. The exact nature of these stresses and their effect on the me- 


chanical properties of the composite has been the subject of numerous investigations. 


Examples of continuum mechanics applied to a thick-walled cylinder model under 
internal pressure or containing a continuous fiber abound in the literature. DeSilva and 
Chadwick [Ref. 1] used a continuous fiber model to investigate dimensional changes due 
to temperature effects in iron-boron and copper-tungsten composites. Their results 1n- 
dicated that matrix vielding occurred first at the fiber interface then spread radially out- 
ward. Dvorak er al. [Ref. 2] used a linear elastic continuous fiber model to determine the 
initial yield surface in the matrix due to temperature changes in the boron-aluminum 
system. They determined that yielding occurred first at the interface and that a tem- 
perature change as small as 10 - 20 °C was sufficient to result in localized vielding in the 
matrix. Hoffman [Ref. 3] used a continuous fiber model to determine thermal stresses 
induced during cooldown in the tungsten fiber reinforced 80N1-20Cr matrix. It was 
conjectured that the magnitude of these stresses would be sufficient to result in plastic 
flow in one or both of the constituents. Garmong [Ref. 4] assumed one-dimensional 
deformation while neglecting shear and poisson effects to calculate matrix stress, strain 
and plastic work as a function of temperature change. He reported the residual stresses 
in the matrix as tensile in nature with plastic strains reaching 0.32%. LEarmme et al. 
[Refs. 5,6] used a continuum mechanics based approach to investigate the stress states 
around spherical, misfitting particles of high symmetry under the assumptions of a von 
Mises vield criterion, linear and power-law strain hardening, and incremental plasticity. 
Arsenault and Fisher [Ref. 7] used the results of electron microscopy and electron 
diffraction experiments to conclude that thermal residual stresses in the matrix were 
caused by CTE mismatch and that matrix strengthening was due to an enhanced dislo- 
cation density in the matrix immediately adjacent to the fiber. It was postulated that 
thermal residual stresses had little effect on matrix strength. In a companion study, 
Arsenault [Ref. 8] compared experimental results with continuum mechanics predictions 


and determined that the presence ofa discontinuous fiber reinforcement resulted in levels 


of strengthening much greater than predicted bv theorv. The cause of this discrepancv 
was an enhanced dislocation density due to the CTE mismatch between the materials. 
Arsenault and Taya [Refs. 9,10] used X-ray diffraction results and yield strength data in 
tension and compression to determine the extent of thermal residual stresses in a SiC 
whisker reinforced 6061 Aluminum composite. The experimental results were compared 
with analytical results obtained using the modified shear-lag theory and the Eshelbv 
equivalent inclusion method. The results indicated that a tensile residual stress existed 
in the matrix on a volume average basis, although tensile and compressive components 
were distributed throughout. Both theoretical methods underpredicted the magnitude 
of strengthening; however the Eshelby method predicted the correct qualitative trends 
for low fiber aspect ratios. The shear-lag theory was seen as a crude approximation even 
with the modification by Nardone and Prewo [Ref. 11] to account for tensile load 
transfer at the fiber tip. A prismatic punching model was developed by Arsenault and 
Shi [Ref 12] to attempt to predict the increase in dislocation density in the matrix on 
cooldown due to CTE mismatch between the constituents. In an exhaustive study, 
Withers er al, (Ref. 13] used the Eshelby method to provide a rigorous theoretical basis 
for the prediction of mechanical properties in discontinuous fiber composites. The 
Eshelby model was modified to account somewhat for inclusion shape, geometry and 
distribution and a full range of elastic, thermo-elastic and plastic composite properties 
were investigated. The objective of the work was to provide a simple, reliable method, 
based on the Eshelby theory, to analyze internal stresses in two-phase materials such as 
metal matrix composites. The mean fiber and matrix stresses were determined as linear 
functions of applied load and fiber/matrix misfit. The model was used to characterize 
composite properties such as elastic modulus, CTE, matrix-to-fiber load transfer, ther- 


mal residual stress levels and interfacial bond strength. Derby and Walker (Ref. 14) 


attempted to characterize the contribution of thermal strains to matrix strength using a 
combination of electron microscopv and hardness tests. 

The use of numerical solution techniques to determine the effect of residual stresses 
on a microstress level have gained wide appeal among researchers. Dutta et al. (Ref. 
15] used a combination of continuum methods with finite element analvsis to model the 
deformation zone arising from differential thermal contraction. Models were con- 
structed for periodic arrays of spherical and cylindrical reinforcements. In a follow-on 
studv, Dutta [Ref. 16] used an axisymmetric finite element model to investigate the effect 
of residual stresses in SiC/Al 6061 composites. The model assumed rigid and elastic fi- 
bers with elastic-ideally-plastic matrices. The results indicated that residual stresses are 
responsible for the less than expected matrix-to-fiber load transfer observed in discon- 
tinuously reinforced mietal-matrix composites. Povirk et al. [Ref. 17] conducted a 
parametric study of the effects of fiber spacing, volume fraction and fiber aspect ratio 
on the residual stresses in metal-matrix composites. The composite was modeled as a 
periodic array of perfectly aligned cylindrical cells, each composed of a cylindrical fiber 
surrounded by a matrix shell. The analysis accounted for thermo-elastic response in the 
fiber and thermo-elastic-plastic response in the matrix. A Lagrangian-Jaumann Stress 
rate formulation was used to model the constitutive response of the material. A 
quenching process was simulated by imposing a temperature history on the cell surface 
through solution of the heat equation for a cylindrical bar. The results indicated that 
transverse fiber spacing was the most important microstructural parameter affecting the 
distribution of residual stress and magnitude of plastic deformation in the matrix. Fiber 
aspect ratio was seen to have negligible effect on residual stresses since these stresses 
became independent of axial position a short distance from the fiber corner. The overall 
level of plastic deformation in the matrix was primarily dependent on the fiber volume 


fraction. 


B. STRESS-STRAIN BEHAVIOR OF METAL-MATRIX COMPOSITES 

Several in-depth studies have recently been conducted to investigate the dependence 
of composite stress-strain behavior on material properties and the size, shape and dis- 
tribution of the reinforcing phase. These studies have relied almost exclusively on finite 
element methods to determine the magnitudes of local field quantities and overall de- 
formation behavior. The flexibility afforded by this method enables a wide range of ef- 
fects to be included in the numerical models that cannot otherwise be accounted for 
effectively using traditional continuum mechanics theories. Factors such as non- 
proportional loading and stress concentration zones due to the presence of sharp corners 
in the reinforcement often render these methods intractable. 

Christman er al. [Ref. 18] conducted an experimental and numerical study of defor- 
mation in metal-ceramic composites with the objective of investigating the dependence 
of tensile properties on reinforcement shape, aspect ratio and distribution. The numer- 
ical analysis was conducted using a Lagrangian Kinematic formulation with Jaumann 
Stress rate and strain rate tensors used in the constitutive response. The matrix material 
was modeled as an elastic viscoplastic solid subject to isotropic work hardening. For the 
fibers, elastic and perfectly rigid models were used in the analysis. The material model 
for the composite consisted of a regular, periodic array of perfectly aligned fibers within 
cylindrical unit cells. The numerical study was carried out within the context of an 
axisymmetric analysis to account for fiber aspect ratio, fiber spacing, and volume frac- 
tion effects; a two dimensional plane strain analysis was used to account for variations 
in fiber distribution. Two different sets of boundary conditions were applied to the 
model to simulate geometric constraint effects. The numerical results indicated that 
significant hydrostatic stress states (up to 6 times the composite offset yield strength) 
were developed in the matrix during deformation. Plastic deformation occurred first at 


the fiber tip then spread rapidly throughout the matrix, with the magnitude of the ef- 


fective plastic strain greatlv exceeding the far-field axial strain. The results of the 
parametric study were in qualitative agreement with experimental results. Fiber volume 
fraction was the single greatest factor affecting composite strength. The effects of fiber 
aspect ratio and spacing were strongly influenced by the geometric constraint imposed 
on the model and the subsequent hydrostatic stress states developed. To compare the 
numerical models with experimental tensile data, average values for fiber and cell aspect 
ratio were chosen from optical microscopy results. The model with imposed sidewall 
constraint predicted values of elastic modulus, vield strength and tensile strength far in 
excess of the experimental data. The model without imposed sidewall constraint over- 
predicted the values of elastic modulus and yield strength, but greatly underpredicted the 
tensile strength and work-hardening rate of the actual composite. 

In two companion studies, Tvergaard [Refs. 19, 20 ] applied the finite element 
scheme used bv Christman et al. [Ref. 18] to investigate the parametric effects of volume 
fraction, fiber spacing and fiber aspect ratio Using a unit cell model in which the fibers 
were shifted relative to one another. The purpose of this modeling constraint was to 
account for matrix shearing between adjacent fiber ends. The author postulated that 
this shearing was a main deformation mechanism in the actual composite. Additionally, 
fiber debonding was accounted for using a cohesive zone model that described normal 
and tangential separation as Well as fiber pull-out. The qualitative results for the volume 
fraction, fiber aspect ratio and fiber spacing effects were identical to the previous study, 
however the effect of matrix shearing due to the shifted fiber arrangement lowered the 
magnitudes of the strength predictions into much closer agreement with experimental 
results. The results also indicated that the assumption of rigid fibers was a good ap- 
proximation for modeling discontinuous fiber composites. The results for the fiber de- 
bonding model revealed that interfacial separation occurred first at the fiber corner 


followed by separation along the top and sidewall fiber surfaces, respectively. The au- 


thor predicted that thermal residual stresses would tend to delav the onset of debonding 
as well as increase Coulombic friction effects once fiber pull-out occurred. 

Using three-dimensional finite element models, Levy and Papazian [Ref. 21] investi- 
gated the tensile properties of aluminum matrix composites containing whisker and 
particulate silicon carbide. The SIC fibers were modeled as longitudinally aligned cylin- 
ders in a three-dimensional array. The fibers were assumed transversely aligned in one 
model and staggered in the other. The matrix material was modeled as elastic-plastic 
and the SiC as elastic. The variables selected for study were fiber voluine fraction, fiber 
aspect ratio, fiber spacing and orientation of the fiber to the applied load. The analytical 
results accurately predicted the increase in strain hardening rate and elastic modulus with 
increasing volume fraction and fiber aspect ratio. The models also predicted the exper- 
imentally observed initial decrease and subsequent increase in proportional limit with 
increasing volume fraction. The results of orienting the fiber orthogonal to the applied 
load indicated little dependence on fiber aspect ratio. The nucrostress data were used 
to predict the location and growth of the plastic zone during deformation. The use of 
three-dimensional analysis allowed the authors to account for the presence of adjacent 
fibers more accurately than in the typical axisymmetric model. The initial yielding lo- 
cation in longitudinal loading was sensitive to fiber aspect ratio and fiber spacing. For 
fiber aspect ratios of one, initial yielding in the matrix occurred at the midpoint between 
adjacent fiber end surfaces. For aspect ratios of four and larger, initial vielding occurred 
along the fiber end surface. In transverse loading, initial yielding always occurred be- 
tween adjacent fibers next to the lateral (cylindrical) fiber surfaces. 

Brockenbrough and Suresh [Ref. 22] used three-dimensional cylindrical and plane 
strain models to investigate the sensitivity of plastic deformation behavior of continuous 
fiber composites to variations in fiber cross-section and distribution. The composite was 


modeled as a periodic array of three-dimensional unit cells with the fibers arranged in 


“square” edge-packing and “hexagonal” edge-packing configurations. Fiber distribution 
effects were studied using a plane strain model of 26 cylinders of aspect ratio 1. The 
analytical results indicated that the plastic deformation characteristics of the composite 
in transverse loading were influenced both by the fiber cross-sectional geometry and 
distribution. However, fiber distribution had a more pronounced effect on transverse 
tensile properties. In longitudinal loading, the evolution of matrix triaxiality and 


plasticity overshadowed both fiber geometry and distribution effects. 


C. RESEARCH OBJECTIVES AND OVERVIEW 

While the previously discussed research efforts have shed light on the extent and 
effect of thermal residual stresses on some strengthening mechanisms in metal matrix 
composites, as Well as the sensitivity of composite tensile properties on processing vari- 
ables, a comprehensive study into the effect of these stresses on the overall stress-strain 
response is not available in the current literature. Accordingly, the purpose of this re- 
search is twofold. The first objective is to investigate the effect of the thermal residual 
stresses generated in the matrix during cooldown from the fabrication’solutionizing 
temperature on the uniaxial stress-strain response of the composite. Since these residual 
stresses are expected to be strongly dependent upon fiber volume fraction, fiber aspect 
ratio and fiber spacing, the second objective is to conduct a parametric study of the ef- 
fect of the above factors on composite tensile behavior. The parametric effects of fiber 
volume fraction, fiber aspect ratio and fiber spacing will be considered within the context 
of two-dimensional, axisymmetric finite element analysis and compared for the cases of 
materials with and without residual stresses, respectively. The model system chosen for 
study is a discontinuous SiC fiber reinforced Al 6061-16 matrix composite. 


Chapter II of this thesis details the assumptions and methods used in the numerical 


modeling of the composite as well as the material testing procedures. Chapter III con- 


tains the results of the investigation and compares the findings with current literature. 


Chapter IV summarizes the relevant conclusions of the study. 


li. EXPERIMENTAL APPROACH 


A. NUMERICAL MODELING 

In order to accurately describe the behavior of a composite material numericallv, the 
model must contain an appropriate level of realism in the choice of geometric and ma- 
terial variables. Figure l is an optical micrograph of a 30 vol.% SiC, 'Al 6061 compos- 
ite. The photograph reveals several characteristics of the actual composite that need to 
be considered in attempting to model the material numericallv. First, there is a consid- 
erable degree of randomness in fiber alignment, with some fibers actuallv oriented 
orthogonal to the tensile axis. Takao et al. (Ref. 23) investigated the effects of fiber 
misorientation on the effective elastic modulus of short fiber composites. Using a 
continuum approach based on the Eshelbv method, the authors found that fiber align- 
ment had a greater impact on composite elastic modulus than fiber volume fraction. A 
second modeling consideration is the non-uniformitv in fiber distribution throughout the 
material. This leads to regions in the matrix in which the fibers are clustered together 
resulting in a local volume fraction in excess of the bulk volume fraction of the material. 
Similarly, there are low-fiber concentration regions in which the local mechanical prop- 
erties of the composite are approximately those of the unreinforced matrix. Previous 
investigations [Refs. 18,22] have described the effect of fiber clustering within the context 
of plane strain finite element models. Their results indicated that clustering greatly re- 
duced the resistance to plastic flow in the composite by decreasing the magnitude of 
matrix triaxiality; however, the elastic properties of the composite were not significantly 
affected by variations in fiber distribution. A third modeling consideration is the wide 
range of fiber aspect ratios encountered in the material. A majority of previous studies 


have attempted to account for this effect by using an average fiber aspect ratio in the 
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numerical calculations. An investigation by Takao and Taya [Ref. 24] revealed that 
volume averaging resulted in larger values of aspect ratio than a simple frequency- 
averaging approach. Experimental data was compared with analvtical results based on 
the Eshelby method. For short fiber composites, it was determined that the use of a 
volume-averaged fiber aspect ratio gave predictions of elastic modulus within 7% and 
CTE within 10% of experimental values. Additionally, the authors indicated that the 
wider the range of aspect ratio distribution in the material, the greater the error in the 
use of an average aspect ratio in analytical results. 
1. Finite Element Modeling 

In using the finite-element method to predict macro-mechanical composite be- 
havior, an assumption implicit in the modeling is that the behavior of a microscopic 
portion of the material can adequately represent the behavior of the bulk material. To 
this end, a multipurpose finite-element code, ADINA, was used to predict macrscopic 
mechanical behavior in the form of uniaxial stress-strain curves, and microscopic be- 
havior in the form of von Mises effective stress contours plots. ADINA permitted great 
flexibility in modeling through the use of geometric and materially non-linear solution 
capabilities. Two-dimensional axisymmetric analysis was used in conjunction with an 
Updated-Lagrangian-Jaumann (ULJ) kinematic formulation. This numerical scheme 
uses the principle of virtual work to calculate incremental displacements, which through 
appropriate constitutive relations yield elemental stresses. The non-linear response of 
the material is approximated through an incremental approach that predicts the config- 
uration of the body at t + At with the present configuration at t. In this case, t denotes 
an iterative counter and not time. This scheme, based on the Newton-Raphson 
equations [Ref. 25], applies incremental displacements to the body until the applied ex- 
ternal load is balanced by the sum of the internal nodal forces. The numerical models 


were constructed using two-dimensional isoparametric elements. Second-order Gaussian 
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integration was used to provide suitable numerical accuracv in the solution. The alu- 
minum matrix was characterized bv a bilinear stress-strain curve, the slopes of which 
represent the elastic modulus, E, and strain-hardening modulus, £,, respectively. The 
aluminum matrix was assumed to undergo isotropic work hardening and obey the von 
Mises yield criterion. All material properties were assumed to be constant witlun the 
temperature range considered in the study. This is a highly idealized assumption in the 
case of aluminum; however, for the purposes of this study, room temperature properties 
were used in view of the modeling of the quench process as an instantaneous event. For 


aluminum 6061 in the T6 condition, the relevant thermo-mechanical properties are E 


oa G50MMPa, c, = 276 MPa, v = 0.33 and CTE = 25.0 u KT); the SiC 
fiber values are E = 450 GPa, v = 0.2 and CTE = 4.3 yA”. 

In the present work, two unit cell configurations similar in nature to models of 
previous researchers [Refs. 18, 19 ] are considered to account for the cases of aligned fi- 
ber (Figure 2a) and staggered fiber (Figure 2b) arrangements, respectively. In ene se, 
the material is assumed to be composed of a periodic array of cvlindrical unit cells per- 
fectly aligned with the tensile axis. As shown in Figure 2c, each cell 1s composed of a 
cylindrical fiber of radius R, and length L, surrounded by a matrix shell of radius R, and 
length £,. In order to simplify computational effort, the three-dimensional problem was 
reduced to two-dimensions by using an axisymmetric modeling approach. The 
axisymmetric approximation assumes that a structure with rotational symmetry about 
an axis in the presence of a rotationally symmetric load can be fully described using a 
two-dimensional model with a unit radian width [Ref. 26]. From Figure 2c, it can be 


seen that the axisymmetric model is an appropriate analog to a three-dimensional array 


of hexagonal unit cells [Ref. 19]. 
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Figure 2. Material Representations: Diagram showing the (a) aligned fiber and 
(b) staggered fiber arrangement with associated unit cell configurations 
and (c) the depiction of the fiber embedded in a matrix shell. 
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2. Boundarv Conditions 

In ADINA, displacement-based boundarv conditions are specified in addition 
to a master set of translational and rotational degrees of freedom for the entire finite- 
element model. In axisymmetric analysis, ADINA requires the model to be contained 
in the positive quadrant of the y-z plane. The master degrees of freedom in this case 
allow y and z translation only. Subsequent specifications on displacement boundary 
conditions can modify, but not delete the master degrees of freedom. 

A cylindrical coordinate system was used with displacements u, and u, corre- 
sponding to the z and r axes, respectively. For the aligned fiber model, two sets of 
boundary conditions were employed to simulate the effect of neighboring unit cells 
through the use of geometric constraints. One model was used to simulate a material 
allowed to undergo unconstrained plastic flow. The only boundary conditions required 
for this case are u, = 0 along z = O and u, = O along r = 0. These conditions describe 
a cell in which deformation response is unaffected by the presence of neighboring cells. 
While such a condition is physically impossible, its value as a numerical model is in 
Serving as a comparison to assess the role of varying degrees of geometric constraint on 
composite constitutive response. 

For the aligned fiber model “constrained” by adjacent unit cells, the boundary 
conditions, in addition to those for the “unconstrained” case, require that the displace- 
ments u, on z = £, and u, onr = R, be constrained to preserve the right circular cylin- 
drical shape of the unit cell throughout the deformation history of the material. In this 
case, the cell “senses” the presence of neighboring cells by application of additional 
constraints at the top and lateral surfaces of the cell such that they remain planar 
throughout the quench and tensile loading process. 

The staggered fiber model (Figure 2b) depicts the composite as a periodic array 


of fibers perfectly aligned with the tensile axis, yet shifted relative to one another along 
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the z-direction. This model is complicated by the necessity of including the combined 
effects of two adjacent unit cells in the boundary conditions. Due to the unique ar- 
rangement of the cells, the top surface of one cell 1s constrained to move not only ina 
plane, but also only as much as the bottom surface of the adjacent cell. Since this bot- 
tom surface contains the fiber, this portion of the cell will contract far less during the 
quench. The requirement for top and bottom surface displacements to be equal and 
opposite imposes greater constraint on the cell than in the aligned “with constraint” case. 
Along the cell centerline, symmetry conditions require that u, = 0 alongr = 0. On the 
lateral surface of the cell at r = R,, compatibility requirements in the axial direction re- 
quire that the u, displacements above and below point C be equal and opposite since this 
point is the geometric center of the two cell svstem. In the radial direction, the situation 
is complicated bv the necessitv to constrain the movement of the nodes above and below 
point C asa Me. of the radial displacement of point C (u,) to ensure that the total 
cross-sectional area of the two cell system remain independent of the axial coordinate. 


Tvergaard [Ref. 19] applied this boundary condition by requiring 
(um) + Re}? + {u,(5) + Ry} = Au, + Ro (1) 


where the coordinates y and ó correspond to the local coordinate systems of the top and 
bottom cell surfaces, respectively. The displacement-based boundary modeling capabil- 
ities of ADINA preclude meeting this requirement explicitly, hence radial compatibility 
is relaxed somewhat in the staggered model. 

The staggered model boundary conditions are an important counterpart to the 
aligned models in that the effects of matrix shearing due to overlapping fibers on the 
plastic deformation behavior of the composite can be studied. Due to the non- 
uniformity of fiber distribution in the material, actual composite response will be a 


combination of both aligned and staggered fiber configurations. 


Perfect bonding is assumed along the fiber-matrix interface. This is accom- 
plished in ADINA by specifying that the fiber and matrix elements adjacent to the 
interface share the boundary nodes, thereby precluding the possibility of separation 
during the loading sequence. 

3. Applied Loading 

The loading scheme used in the numerical modeling combines thermal and dis- 
placement loads in a sequence designed to approximate those experienced by the com- 
posite during fabrication and subsequent tensile testing. 

The cooldown from the fabrication/solutionizing temperature to room temper- 
ature is accomplished by subjecting all nodes in the unit cell to a 505 °C temperature 
drop in one solution time step. In ADINA this is accomplished by selecting the incre- 
mental time step size to be equal to the solution time period, therebv bringing the nodal 
temperatures from the solutionizing temperature (330 *C) to room temperature (25 *C) 
in one step. The rationale for adopting this approach was to: (1) eliminate the effect of 
quench rate from the solution and (2) eliminate size and bulk thermal property effects 
necessary when attempting to impose a thermal history at the surfaces of the unit cell. 
While this analysis does not take into account factors such as thermal softening due to 
heat generation during plastic deformation, it was assumed that the dominant effect on 
the composite from the fabrication cycle is the generation of an enhanced dislocation 
density in the vicinity of the fiber-matrix interface that is primarily a factor of CTE 
mismatch. The resulting nodal displacements from the quench were stored in a restart 
file to be used as initial nodal displacements during tensile loading analysis. The use of 
restart analysis allows great flexibility in modeling by permitting the use of a variety of 
loading conditions in a sequence, thereby providing a more realistic simulation of real 


processes. 
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Tensile loading of the model was accomplished bv applving an axial displace- 


ment of 
AL, = (e° — I)L, (2) 


where e is the applied true strain. For the aligned model, AL, was applied to the top 





. : | AL 
surface of the unit cell; in the staggered case, a displacement of 5 
top and bottom surfaces, respectively. The composite stress-strain curves Were con- 


— was applied to the 


structed by applying a AL, corresponding to a selected value of e. The average axial 


stress was calculated using a volume-average approach as follows: 


=p [dav + a | otal, (3) 
C C 
Where superscripts c, f and m correspond to composite, fiber and matrix values, respec- 
tively. In this equation, each c,, represents an average of the four stress values obtained 
from the Gaussian integration points within each finite-element. The average elemental 
Stresses are then multiplied by their corresponding elemental volumes and summed over 
the entire volume of the composite. 
4. Grid Independence 

In order to obtain accurate solutions from a numerical model, convergence crl- 
teria for the linearization process and the mesh geometry must be established. In 
ADINA, the user has the option of selecting the iterative scheme for linearization, the 
basis for convergence (displacement, energy, etc.) and the value for the convergence 
criterion. For the Full-Newton iterative scheme selected for use in this study, out-of- 


balance energy was the basis for convergence and the required tolerance was set to 


10-6. The mesh for each model was uniformly refined until 2 successive values of o%, at 
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e= Q (the residual axial stress) gave results that agreed to within 10%. The residual axial 
Stress value was found to be the factor most sensitive to changes in mesh size and, 
therefore, was determined to be the most suitable parameter to establish grid independ- 
ence. In most cases, the numerical models achieved grid independence with approxi- 
mately 250 elements. 

3. Model Implementation 

The numerical models discussed above were modified for the particular variable 
whose effect was being evaluated. The models used to approximate the 30 vol.% ex- 
perimental data were assigned values of fiber and cell aspect ratios of 5 and 6, respec- 
tively while the 20 vol.% models used values of 3 for both fiber and cell aspect ratios. 
For volume fraction effects, the fiber and cell aspect ratios were both fixed at 3 and nu- 
merical results were carried out for volume fractions of 1090, 2090 and 3090. For the 
fiber aspect ratio study, the fiber volume fraction was fixed at 30%; the fiber and cell 
aspect ratios were set equal to one another to eliminate spacing effects and numerical 
results carried out for values of 1, 3, 6 and 10. For the fiber spacing study, the fiber 
volume fraction was again fixed at 30%; the cell aspect ratio was fixed at 6 and numer- 
ical results were carried out for fiber aspect ratios of 3, 6 and 10. 

Contours of von Mises effective stress were generated to determine the location 
of initial yielding in the matrix and the subsequent growth of the plastic zone due to 
loading. These models were developed using approximately twice the number of ele- 
ments used in the constitutive response models to ensure clarity and resolution of the 


contour curves. 


B. MATERIAL TESTING 
The material testing portion of this study primarily focused on the tensile testing of 
both unreinforced and fiber-reinforced Al 6061 to extract mechanical properties for use 


in the numerical models , and in the case of the fiber-reinforced aluminum, to serve as 
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a comparison to assess the accuracy of the numerical predictions of stress-strain re- 
sponse. The composite material was procured from Advanced Composite Materials 
Corporation (ACMC), Greer, SC, and was produced using powder metallurgy (P/M) 
techniques. Extruded billets were subsequently produced from the P/M compacts for 
20 vol.% SiC/Al 6061 at a reduction ratio of 11:1 and 30 vol.% SiC/Al 6061 at a re- 
duction ratio of 6.5:1. The composite matrix was heat treated to the T6 condition by the 
manufacturer. Tensile test specimens with 1.0” gage lengths were cut from the extruded 
stock such that the longitudinal axis of the sample was aligned with the extrusion di- 
rection. Unreinforced Al 6061-T6 tensile specimens were machined from plate stock. 

The tensile specimens were mounted in a screw-driven Instron machine and pulled 
at rates of 1.0 mm/min for the aluminum samples and 0.127 mm, min for the composite 
samples. Clip-on extensometers were used to record engineering strain values during the 
test. The curves for the Al 6061 samples were used to determine values of yield strength 
and strain-hardening modulus for input into the numerical models. The value of the 
strain-hardening modulus was determined by fitting a curve from the offset yield point 
to the stress-strain curve within the expected strain range of the composite, in this case 
5%. A more accurate value would be obtained by using an unreinforced powder sample 
from the same lot that the composite was produced, followed by heat treatment together 
with the composite samples to ensure uniformity in material processing. 

An optical micrograph was obtained from a sample of the 30 vol.% extruded stock 
to determine average values of fiber aspect ratio and cell aspect ratio as well as to serve 


as a microstructural comparison When developing the unit cell models. 
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HI. RESULTS AND DISCUSSION 


A. COMPARISON OF ALIGNED AND STAGGERED FIBER MODELS 

The sensitivity of composite constitutive response to imposed modeling constraint 
can be seen in Figure 3. As shown in Figure 3a, the unconstrained aligned-fiber model 
for a composite with 30 vol.% SiC whiskers of aspect ratio 5 predicts a Young's modulus 
(E) of 148 GPa, a proportional limit (PL) of 181 MPa and a yield strength (YS) of 276 
MPa. The aligned-fiber model with constraint, in comparison, predicts the composite 
E, PL and YS values as 163 GPa, 292 MPa and 429 MPa, respectively. The constrained 
model also predicts a substantially higher work-hardening rate than the unconstrained 
model. Figure 3a also shows the experimental data for a 30 vol.% metal-matrix com- 
posite (MMC). By comparing analytical results with experimental data, it is evident that 
both the constrained and unconstrained models overestimate composite stiffness slightly. 
This can be attributed to the fact that the models assume perfectly aligned, parallel fi- 
bers. From Figure I, the fibers in the experimental material are clearly neither aligned 
nor parallel. Additionally, some deviation is expected due to the mismatch between the 
fiber aspect ratio and fiber spacing assumed in the model and that occurring in the actual 
composite. The unconstrained model was found to underestimate the composite yield 
strength and work-hardening rate considerably. This is expected since the matrix in the 
unconstrained model is unable to “sense” the presence of adjacent unit cells and, there- 
fore, is allowed to undergo relatively unconstrained plastic flow. The constrained model, 
on the other hand, is a much more realistic representation of an actual composite, al- 
though it predicts an appreciably higher yield strength and work-hardening rate. This, 
as discussed earlier, is attributable to the considerable deviation of the actual MMC 


from the idealized fiber configuration assumed in the model. 
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Figure 3. Modeling Constraint Effects: Comparison of aligned and staggered 
models with a (a) 30 vol.% tensile specimen and (b) 20 vol.% tensile 
specimen, respectively. 


me 


As seen in Figure 3a, the staggered model provides a closer approximation of the 
experimental data, despite the partially relaxed radial compatibility conditions (section 
11.A.2). The model predicts stiffness, proportional limit and work-hardening rate values 
which are reasonably close to the actual composite, but overestimates the yield strength. 

The above discussion suggests that the degree of fiber alignment along the loading 
axis and the extent to which the adjacent fibers are shifted with respect to each other 
has as large effect on the accuracy of the analytical predictions. The degree of fiber 
alignment increases with increasing extrusion ratio, and hence it is expected that the 
agreement between analytical results and experimental data will improve with larger 
extrusion ratios. In the limiting case where the extrusion ratio is very high and all the 
fibers are perfectly aligned along the loading axis, the predictions of the aligned model 
are expected to match the experimental data reasonably well, although some over- 
prediction is likely since sone degree of lateral misalignment is probable. In that case, 
the staggered fiber model is probably a more reasonable approximation, provided the 
radial boundary conditions of the model are adjusted appropriately. With the present 
radial boundary conditions, however, the staggered model is expected to somewhat 
underestimate the mechanical properties (PL,YS) of a composite with a very large 
extrusion ratio (i.e. a high degree of fiber axial alignment). 

Figure 3b rates the effect of the imposed modeling constraints compared with 
20 vol.% experimental data. Both aligned and staggered models give reasonable pred- 
ictions of stiffness, with experimental values of proportional limit, vield strength and 
work-hardening rate lying between the values predicted by the numerical models. As 
expected, because of the higher extrusion ratio of the 20 vol.% sample (11:1, compared 
to 6.5:1 for the 30 vol.% sample), the predictions of the constrained aligned-fiber model 
are in closer agreement with the experimental data than at 30 vol.%. The staggered 


model now underpredicts the composite response somewhat due to the relaxed radial 
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compatibility condition. From the above, it is reasonable to assume that the aligned- 
fiber model is a fair approximation of a real composite provided the extrusion ratio is 
very large. The staggered model with the present boundary conditions, however, while 
yielding reasonable predictions for composites with low to medium extrusion ratios, is 
not a correct representation of the constraint mechanisms operative in an actual MMC. 
Hence, the aligned-fiber model was the model of choice for the parametric study pre- 
sented in subsequent sections. 

The close agreement between analytical results and experimental data at low fiber 
volume fractions is well represented in the literature [Refs. 18,19]. The majority of pre- 
vious studies have made numerical comparisons with control systems of volume frac- 
tions from 5% to 20%. The greater level of error at higher volume fractions can be 
attributed to the modeling process; at high volume fractions, the reduced lateral and 
longitudinal spacing between the fiber and cell result in increased levels of matrix con- 
Straint, and hence have a significant effect on the plastic flow behavior of the material. 
Hence a slight error in the estimation of fiber aspect ratio and fiber spacing translates 
into a large error in the prediction of composite response. For low volume fraction 
models, the boundaries of the unit cell are relatively free of these constraint effects and 
the numerical results more closely represent experimental results. 

Another factor affecting the trend of the numerical results is the difference in mean 
fiber aspect ratio between the 20 and 30 vol.% samples. The 20 vol.% sample was 
processed at a reduction ratio of 11:1, whereas the 30 vol.% sample was processed at a 
ratio Of 6.5:1. The end result of this difference is that the 20 vol.% sample experiences 
more fiber breakage during processing and subsequently exhibits a lower value of aver- 
age fiber aspect ratio (4;= 3) than the 30 V% sample (4,= 5). Christman et al. [Ref. 
18] observed that the use of particulate-reinforced models (A,=1) produced better 


agreement with experimental data than using whisker-reinforced systems. The reasons 
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for increased accuracy of the numerical models were (1) particulate misalignment with 
the tensile axis had a negligible effect on composite response for 4,=1, (2) the 
particulate morphology produced less plastic flow constraint in the composite with the 
result that the idealized cell model more closely emulated the behavior of the actual 
material and (3) clustering effects in particulate-reinforced composites had less of an ef- 


fect on the stress-strain response than in whisker-reinforced materials. 


B. EFFECT OF FIBER VOLUME FRACTION 

Figure 4 illustrates the effect of increasing volume fraction on the constitutive re- 
sponse of the composite. The numerical models predict the experimentally observed 
trends of increasing values of stiffness, vield strength and work-hardening rate with in- 
creasing volume fraction. The primary mechanisms affecting composite behavior in this 
case are increased plastic flow constraint, increased dislocation density and enhanced 
matrix-to-fiber load transfer with increasing reinforcement volume fraction. Christman 
et al. (Ref. 18] observed that as volume fraction increased, the magnitude of the 
hydrostatic stresses in the matrix dramatically increased producing greater resistance to 
plastic flow, and hence higher work-hardening rates. Another trend predicted by the 
models is the experimentally observed behavior of decreasing and subsequent increasing 
values of proportional limit with increasing volume fraction. Levy and Papazian (Ref. 
21] reported this effect as a function of localized stress concentrations in the matrix at 


low levels of reinforcement addition. 


C. EFFECT OF FIBER ASPECT RATIO 

The dependence of composite constitutive response on fiber aspect ratio can be seen 
in Figure 5. The numerical results predict an increase in stiffness, proportional limit, 
yield strength and work-hardening rate with increasing aspect ratio. Additionally, the 
curves indicate that the strengthening effect becomes less pronounced as the fiber aspect 


ratio increases. The increased strength of composites with larger aspect ratios can be 
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Figure 4. Volume Fraction Effects: (a) Curves of predicted composite response 
for reinforcement volume fractions of 0%, 10%, 20% and 30%, respec- 
tivelv and (b) schematic of model configurations. 
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Figure 5. Fiber Aspect Ratio Effects: Curves of composite response for fiber as- 
pect ratios of 1,3,6 and 10, respectively. 
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attributed to increased load transfer between the matrix and the longer fibers. Piggott 


[Ref. 27] described this strengthening effect in terms of a critical aspect ratio (1,,,,), below 


ent 
which the load transfer between the matrix and fiber fail to produce the maximum fiber 
axial stress possible. Figure 6 illustrates this effect. For fiber aspect ratios above the 
critical value, surface shear transfer from the matrix can result in maximum load transfer 
to the fiber. For aspect ratios equal to the critical value, portions of the fiber can still 
experience maximum load transfer, although the average stress in the fiber 1s lower. 
When the aspect ratio 1s reduced to below the critical value, the load transfer is such that 
the resulting values of axial stress and average stress in the fiber are far below the pre- 
vious cases. In this case, the relative inability of the matrix to transfer load to the fiber 
results in the matrix bearing a substantially greater portion of the applied load, leading 
to lower values of composite stiffness, yield strength and work-hardening rate. It 1s 
generally believed that matrx-to-fiber load transfer is not a major strengthening mech- 
anism in discontinuously-reinforced MMC's [Refs. 28 ,29], which are thought to be 
strengthened primarily by increased matrix dislocation density. However, as shown in 
Figure 7, the average axial stress supported by the fiber increases with both increasing 
fiber volume fraction and fiber aspect ratio. It is evident that load transfer from the 
matrix to the fiber is an important factor in the strengthening of fiber-reinforced 
MMC's. This is further supported by the increase in Young's modulus observed both 
analytically and experimentally with increasing fiber volume fraction. If an increased 
dislocation density were the only strengthening mechanism, no increase in stiffness 
would be observed on increasing the reinforcement content of the composite. 

The reduced degree of additional strengthening with increasing aspect ratio at higher 
values of fiber aspect ratio is observed only for aligned fibers. The trends observed in 
other studies [Refs. 19,21] using the staggered-fiber arrangement indicate that the level 


of strengthening in the composite increases with increasing aspect ratio. This is contrary 
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to the trend observed in the aligned-fiber model. One explanation of the difference in 
behavior is due to the effect of fiber overlap in the staggered model. As fiber aspect ratio 
increases, the degree of fiber overlap increases producing greater geometric constraint 


and therefore increased resistance to plastic flow. 


D. EFFECT OF FIBER SPACING 

The sensitivity of composite constitutive response to variations in fiber spacing is 
shown in Figure 8a. The results indicate that changes in fiber spacing have little effect 
on composite elastic properties; however, the resistance to plastic flow in the material is 
significantly affected bv close end-to-end and side-to-side fiber spacing. The curve re- 
presenting close end-to-end spacing (4,=3) exhibits the greatest increase in work- 
hardening rate due to the relatively small amount of matrix material separating adjacent 
fiber ends. This results in a substantial increase in load transfer across the fiber end 
surface, and hence increased resistance to plastic flow. The curve representing close 
side-to-side fiber spacing (4,= 10) exhibits a somewhat lower work-hardening rate than 
observed in the close end-to-end spacing case. The lateral constraint imposed on the 
composite in this case, while not as great an effect as for A,=3, still increases the re- 
sistance to plastic flow in the material above that of the composite model without 


spacing constraints (4,= 6). 


E. EFFECT OF THERMAL RESIDUAL STRESSES 

The primarv aim of this studv has been to assess the role of thermal residual stresses 
on the stress-strain behavior of metal-matrix composites. The results contained in this 
section describe this effect in terms of fiber volume fraction and geometric constraint due 
to the modeling of fiber arrangement. In addition, load-bearing characteristics of the 
fiber during deformation is studied in an effort to more fully explain the stress-strain 


behavior of the numerical models. 
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Figure 8. Fiber Spacing Effects: (a) Numerical results for fixed volume fraction 
and fiber aspect ratio with variable fiber spacing and (b) schematic of 
fiber spacing models. 
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Residual Stress Effects 
For Af=3 Ac-3 
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Figure 9. Residual Stress and Volume Fraction Effects: Effects of residual stress 
as a function of (a) volume fraction. Also shown are the stresses sup- 
ported by the fiber as å function of strain for (b) 10 vol.% and (c) 30 
vol.%. 
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dual stress predicts a higher composite stiffness and vield strength, while the model with 
residual stress exhibits a higher work-hardening rate. The results for 30 vol.% indicate 
that both curves exhibit similar behavior in the elastic regime, however the curve with 
residual stress displays a slightly higher work-hardening rate. To explain these results, 
plots of fiber stress-versus-strain were compared with the corresponding composite 
curves. The results indicate that for 10 vol.% the presence of residual stresses inhibits 
the load transfer to the fiber resulting in an overall decrease in composite strength in 
both the elastic and plastic regimes. The results for 30 vol.% indicate the opposite trend. 
The presence of residual stresses enhances the transfer of load to the fiber enabling the 
composite as a whole to exhibit greater strength properties. From Figure 9, it is also 
evident that the rate of load transfer to the fiber with increasing strain is much higher 
in the 30 vol.% composite. While for the 10 vol.% composite, the rate of load transfer 
with and without residual stresses is fairly close, the 30 vol.% exhibits a much higher rate 
initially in the presence of residual stresses, leading to the observed cross-over of the fi- 
ber stress-strain curves. This effect indicates that for the 30 vol.% composite, the pres- 
ence of residual stresses produces a much higher rate of load transfer from the matrix 
to the fiber resulting in the observed increase in composite mechanical properties. 
Figure 10 illustrates the role of thermal residual stresses as a function of imposed 
geometric constraint. In this case, the aligned and staggered models are compared ata 
volume fraction of 30% and equal cell and fiber aspect ratios of 3. The results indicate 
that the aligned model with residual stress predicts greater strengthening effects in both 
the elastic and plastic regimes, whereas the staggered model exhibits lower stiffness and 
yield strength and a higher work-hardening rate in the presence of residual stresses. The 
fiber stress-strain results for the staggered model indicate that the transfer of load from 
the matrix to the fiber is inhibited in the presence of residual stresses, similar to the re- 


sults for the aligned model at 10 vol.% (Figure 9b), 
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Figure 10. Residual Stress and Modeling Constraint Effects: Effects of residual 
stress as a function of (a) geometric constraint. Curves of fiber stress- 
strain response for a (b) 30 vol.% aligned and (c) 30 vol.% staggered 
model. 
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F. PLASTIC ZONE INITIATION AND GROWTH 

In the previous section, the role of thermal residual stresses as functions of volume 
fraction and modeling constraint were discussed in terms of their effect on the overall 
stress-strain response of the composite. In this section, the role of residual stresses on 
the initiation and growth of the plastic zone during thermal and subsequent tensile 
loading of the composite 1s studied with reference to the effect of reinforcement volume 
fraction, modeling constraint and fiber spacing on the effective stress state in the com- 
posite in an effort to identify the main deformation characteristics operative in an actual 
composite. 

I. Thermal Residual Stress Effects 

Figure 11 shows a 30 vol.% aligned model at selected stages of deformation. 

The contours of residual effective stress (e =0) reveal that a plastically deformed zone 
(0,> 275 MPa) exists along the entire fiber-matrix interface with the highest stress 
states occurring along the fiber corner. The plot also indicates that the entire area of the 
matrix from the lateral interface to the cell boundary has undergone plastic deformation. 
The stress levels in the fiber indicate a compressive stress zone in excess of 500 MPa 
concentrated along the fiber corner and extending along the end surface. This high 
stress zone quickly dissipates to below 250 MPa a short distance from the end surface 
of the fiber. At £=0.001, the plastic zone along the lateral surface of the interface begins 
to expand in the radial direction. The matrix region adjacent to the end of the fiber ex- 
periences a stress relaxation as load is transferred from the matrix to the fiber across the 
fiber tip and end surface. This effect can be seen in the growth of the region of the fiber 
tip stresses above 500 MPa,indicative of the fiber bearing a greater portion of the applied 
loading at e=0.001. At e=0.0045, the entire fiber is stressed above 500 MPa and the 


plastic zone has further developed in the area between the lateral interface and cell 
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Figure 11. Plastic Zone Development-Residual Stress Effects: Contours of eflec- 
tive stress for a 30 vol.% aligned model with A, = A, = 6. 
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boundarv and has begun to expand to the matrix region above the fiber as evidenced 
bv the 2/5 MPa contour lines. 

The contours for the 30 vol.% model without residual stresses illustrates the 
fundamental difference in behavior of the composite with and without residual stresses. 
In the model without residual stress, plastic deformation begins in the matrix region 
above the fiber and at the fiber corner. At increased load levels, the plastic zone expands 
in the axial direction into the matrix material above the fiber, followed by radial growth 
along the lateral fiber interface out to the cell boundary. Thus, when residual stresses 
are absent, the matrix region next to the lateral (cylindrical) fiber surface is the last to 
deform plastically, suggesting that the matrix in this region 1s able to transfer a large 
portion of the applied load to the fiber. Little, if any, load transfer occurs across the 
fiber ends in this case. When thermal residual stresses are present, on the other hand, 
the fibers are in axial compression prior to loading, and a large degree of load transfer 
takes place across the fiber ends, as evidenced by the initial stress relaxation in the ma- 
trix and the subsequent development ofa small plastic zone upon external loading. The 
increased load transfer across the fiber end surface, however, 1s at the expense of some- 
what reduced load transfer across the lateral portion of the fiber-matrix interface, as 
evidenced by the increase in the effective stress in the plastic zone adjacent to the inter- 
face from the outset. 

2. Volume Fraction Effects 

Figure 12 shows contours of effective stress for a 10 vol.% aligned model. The 
results clearly illustrate the effect of varying reinforcement volume fraction on the de- 
formation characteristics of the composite. The absence of fiber spacing constraints in 
the 10 vol.% numerical model results in a residual stress field in the matrix characterized 
by larger levels of plastic deformation in the vicinity of the fiber-matrix interface than 


for the 30 vol.% model. The fiber exhibits an extremely high compressive residual stress 
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Figure 12. Plastic Zone Development- Volume Fraction Effects: Contours of ef- 
fective for a 10 vol.% aligned model with A,= A,= 6. 
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characterized by a 500 MPa zone in the vicinity of the fiber corner, a low stress region 
in the material below the end surface and an adjacent high stress region (above 500 
MPa) extending to the midpoint of the fiber. This 1s contrasted with the 30 vol.% model 
in which the region below the top surface of the fiber 1s relatively unstressed (below 250 
MPa). Physically, the 10 vol.% model has a larger volume of matrix material sur- 
rounding the fiber which, upon quenching, encounters less restriction in movement than 
in the 30 vol.% model. The matrix, therefore, is free to deform plasticallv to a greater 
extent in the vicinity of the interface. The smaller end surface area of the fiber precludes 
the level of load transfer across this surface exhibited in the 30 vol.% model. The end 
surface of the fiber fails to accommodate the stress build-up in the matrix in this region 
resulting in increased levels of matrix plasticity. The extent of matrix material above the 
fiber in the 10 vol.% model, being greater than in the 30 vol.% model, implies that for 
a given applied strain, the matrix in the 10 vol.% model must transfer the subsequent 
displacements a greater distance prior to encountering the fiber. This results in lower 
levels of load transfer, and hence lower levels of composite strength in the 10 vol.% 
model. 

At low levels of applied strain (e = 0.001), the matrix material adjacent to the top 
surface of the fiber experiences a stress relaxation as the fiber accommodates the applied 
load. During this stage, the plastic zone grows only slightly as the majority of the load 
is transferred to the fiber. At e=0.0045, the plastic zone experiences both radial and 
axial growth with the greatest expansion occurring in the radial direction as in the 30 
vol.% model. The matrix in the region well above the fiber interface exhibits a umform 
stress state of 275 to 280 MPa. Noteworthy is the lack of high and low stress zones 
evident in the matrix region above the fiber in the 30 vol.% model. In the 30 vol.% 


material, these zones presumably result from increased levels of matrix hydrostatic stress 
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due to the higher volume of reinforcement, although this premise has not been investi- 
gated. 

The 10 vol.% model in the absence of residual stresses exhibits the same axial 
growth of the plastic zone as in the 30 vol.% model. The fiber in this case, however, 
accommodates a greater level of load transfer from the matrix than with residual stresses 
present. This provides one explanation of the results presented in Figure 9 (i.e. the de- 
crease in yield strength due to the presence of residual stresses in the 10 vol.% MMC 
versus the increase in vield strength observed in the 30 vol.% MMC in the presence of 
residual stresses). In the 10 vol.% composite, the degree of load transfer across the fiber 
ends in the presence of residual stresses is much less than in the 30 vol.% composite. 
Hence, while the presence of residual stresses vields greater composite strength in the 
30 vol.% model, it reduces the strength in the 10 vol.% composite since the decrease in 
load transfer across the lateral fiber surface is not offset by a corresponding increase in 
load transfer across the fiber end surface. 

3. Effects of Modeling Constraints 

Figure 13 shows contours of effective stress for a 30 vol.% staggered fiber 
model. The premise for studying the staggered fiber arrangement was the expectation 
that significant stress states would be developed in the matrix in the region between 
overlapping fiber ends, resulting in an important deformation mechanism operative in 
an actual composite material. A comparison of Figure ll and Figure 13 reveals a 
profoundly different pattern of plastic zone development. The plot of residual effective 
stress reveals a complicated deformation pattern composed of (1) plastic zones adjacent 
to the fiber interface of larger magnitude than existing in the 30 vol.% aligned model, 
(2) larger levels of plastic deformation in the matrix region between the lateral fiber 
interface and the cell boundary and (3) regions of plastic deformation in the matrix well 


away from the fiber interface that is a result of the boundary conditions imposed on the 
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Figure 13. Plastic Zone Development-Constraint Effects: Contours of effective 
stress for a 30 vol.% staggered model with /t, < 4,< 6. 
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numerical model for the staggered fiber arrangement. A fundamental deformation 
characteristic of this model in the presence of residual stresses is the growth of the plastic 
zone in the matrix region above the interface first before the radial growth common in 
the aligned models occurs. At higher strain levels (e = 0.0045), the plastic regions inter- 
sect followed by gross plasticity in the matrix. 

The contour plots of the staggered model without residual stresses exhibits the 
axial growth of the plastic zone common in the aligned models. The load transfer to the 
fiber occurs more readily in this case than in either the staggered model with residual 
stress or the 30 vol.% aligned model. This factor accounts for the observed increase in 
yield strength and stiffness over the same model with residual stresses present 
(Figure 10a). The observed increase in work-hardening rate in the staggered model with 
MM dual stresses can be attributed to the development of hvdrostatic stress pockets in the 
matrix (Ref. 21] along the cell boundary that offer an increased resistance to plastic flow 
over the model without residual stresses. 

4. Fiber Spacing Effects 

Figure 14 illustrates the effect of fiber spacing in the presence of residual 
stresses on the development and growth of plastic deformation in the composite matrix. 
The results indicate that side-to-side fiber spacing has the greatest impact on the level 
of residual stress in the matrix as shown in Figure 14a. The plastic zone for this model 
extends throughout the entire matrix with the highest stress levels occurring at the 
interface. In contrast, the models with A,= 6 and 10, respectively exhibit a plastic zone 
confined to the region near the fiber interface. Figure 14c shows the stress distributions 
in the composite with close end-to-end spacing. The resulting deformation pattern 
reveals a plastic zone with larger radial extent than the A, < 6 model due to the in- 
creased lateral spacing. The plastic zone on the top surface of the fiber is confined to a 


relatively thin layer of material with an adjacent low-stress zone in the material at the 
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Figure 14. Plastic Zone Development-Fiber Spacing Effects: Contours of effec- 
tive stress for 30 vol.% aligned fiber models with cell aspect ratios of 6 
and fiber aspect ratios of (a) 3, (b) 6'and (c) 10, respectively. 
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relativelv thin laver of material with an adjacent low-stress zone in the inaterial at the 
top of the cell. The stress states in the fiber are a maximum in the ,1,< 3 model followed 
by the 4,76 and 10 models, respectively. These results help explain the stress-strain 
behavior observed in Figure 8. Even though the material with close side-to-side spacing 
initially exhibits greater levels of plastic deformation throughout the matrix, the close 
proximity of adjacent fiber end surfaces produces a more pronounced effect on the 


composite work-hardening rate. 
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IV. CONCLUSIONS 


In this study, the uniaxial stress-strain response of an SiC,/Al 6061 composite was 


investigated to assess the role of thermal residual stresses on the constitutive behavior 


of the material. The following summarize the relevant findings: 


L. 


tv 


Elastic-plastic and thermo-elastic-plastic analysis using finite-element models of 
fiber-reinforced metal-matrix composites can be used to predict uniaxial stress- 
strain behavior given the thermo-mechanical properties of the constituent phases 
as a starting point. 


Aligned and staggered fiber models are both necessary to describe the fundamental 
deformation mechanisms operative in the composite material. 


. The stiffness, yield strength and work-hardening rate of the composite increase with 


increasing fiber volume fraction and fiber aspect ratio. 


Variations in fiber spacing have little effect on composite elastic properties, but 
significantly affect the work-hardening rate for close side-to-side and end-to-end 
spacing. 


The presence of residual stresses significantly impact the stress-strain behavior of 
the composite, principally by affecting the load transfer characteristics between the 
matrix and fiber. 


The presence of residual stresses affect the initiation and continuation of plastic 
deformation in the matrix by producing a plastic zone adjacent to the fiber-matrix 
interface that, in general, expands radially outward from the lateral surface of the 
fiber to the cell boundary as loading increases. The absence of residual stress, in 
contrast, results in the simultaneous development of plastic zones in the matrix 
above the fiber end and along the fiber corner that expand axially as load increases. 


The lateral spacing of adjacent fibers is an important factor affecting the magnitude 
and distribution of residual stresses in the composite. 
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APPENDIX A. SAMPLE ADINA INPUT FILE 


* THIS INPUT FILE SIMULATES THERMAL AND TENSILE LOADING OF A 

* METAL-MATRIX COMPOSITE. THE PROGRAM INCLUDES COMMENTS AT 

* APPROPRIATE LOCATIONS TO ILLUSTRATE SOME OF THE IMPORTANT 

* FEATURES OF ADINA IN NON-LINEAR ANALYSIS. 

de 

% THE FOLLOWING BLOCK OF CODE CONTAINS ADINA CONTROL INFORMATION 
DATABASE CREATE 

CONTROL ORIGIN=UPPERLEFT PLOTUNIT=CM 

FILEUNITS LIST=6 LOG=6 ECHO=6 


* THE FOLLOWING BLOCK OF CODE CONTAINS SPECIFICATIONS FOR MODEL DEGREES OF 
* FREEDOM, SOLUTION TIME PERIOD LENGTH AND TIME-STEP SIZE. THE FIRST 

* "MASTER" COMMAND LINE IS USED IN ELASTOPLASTIC ANALYSIS; THE SECOND 

* COMMAND SIMULATES THE QUENCH PROCESS AND THE THIRD SIMULATES THE 

æ TENSILE LOADING FOLLOWING THE QUENCH. 

ve CAN BE USED IN DURING A GIVEN RUN, THE OTHER MUST BE COMMENTED OUT USING 
a x 

MASTER IDOF=100111 ITP56=0 NSTEP=6 DT=1.0 

MASTER IDOF=100111 ITP56=2 NSTEP=1 DT=1.0 

“MASTER IDOF=100111 ITP56=2 NSTEP=6 DT=1.0 

* THE FOLLOWING BLOCK OF CODE CONTAINS THE KINEMATIC AND ITERATION METHOD 
* SPECIFICATIONS.  AUTOMATIC-ATS IS A FEATURE USED IN NON-LINEAR ANALYSIS 
w TO ASSIST IN ACHIEVING THE REQUIRED TOLERANCE FOR THE OUT-OF-BALANCE 

ve ENERGY NORM. 

KINEMATICS DISP=LARGE STRAINS=LARGE 

ITERATION METHOD=FULL-NEWTON LINE-SEARCH=YES 

AUTOMATIC-ATS N=10 

* 


nA vn 


* THE FOLLOWING COMMANDS RELATE TO THE PROGRAM OUTPUT. COMMENT IN THE 
x "WORKSTATION" COMMAND IF A HARDCOPY OF THE MODEL CONFIGURATION IS 

%* DESIRED. 

HEADING 'VF=30% AF=3 AC=3' 

PRINTOUT VOLUME=MAX 

PORTHOLE FORMATTED=YES FILE=60 

“WORKSTATION SYSTEM=4 DEVICE=0 OPTION=0 

Je 


% THE FOLLOWING CODE SPECIFIES THE NATURE OF THE APPLIED TEMPERATURE LOAD. 
% NUMBERS TO THE LEFT ARE TIME VALUES AND THE NUMBERS TO THE RIGHT ARE THE 
% VALUES OF TEMPERATURE CHANGE FROM THE INITIAL TEMPERATURE OF THE MODEL. 
TIMEFUNCTION 1 

0.0 0.0 

1.0 2505.00 


% THE FOLLOWING TIME FUNCTIONS REPRESENT THE DISPLACEMENT LOADING SEQUENCE 
* THE NUMBERS TO THE RIGHT ARE DISPLACEMENTS FROM THE REFERENCE LENGTH OF 
* THE MODEL THAT CORRESPOND TO A GIVEN LOGARITHMIC STRAIN TO BE USED IN 

* THE STRESS-STRAIN CURVE. NOTE: THE DISPLACEMENTS FOR THE MODEL WITH 

%* RESIDUAL STRESSES START AT TIME=2.0 SINCE THE QUENCH PROCESS OCCURRED 

* FROM TIME 0 to 1.0. IN ADDITION, THE DISPLACEMENTS ARE INITIALLY 

* NEGATIVE DUE TO THE QUENCH PROCESS. THE SUBSEQUENT DISPLACEMENTS 
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* ARE ADDED TO THE NEW REFERENCE. THE DISPLACEMENT SPECIFICATIONS 

* IN ADINA ARE ALWAYS ABSOLUTE DISPLACEMENTS, MEANING THAT A RELATIVE 
* CHANGE IN DISPLACEMENT CANNOT BE INPUT. 

TIMEFUNCTION 2 

* DISPLACEMENTS FOR ALIGNED MODEL WITH RESIDUAL STRESSES 


0. 0 0.0 

250 -. 015496 
340 -. 013248 
4.0 -, 00875 

350 0. 00702 

6.0 0.01124 

7.0 0.06154 


* 


TIMEFUNCTION 3 


* DISPLACEMENTS FOR ALIGNED MODEL WITHOUT RESIDUAL STRESSES 


0.0 0.0 

10 0. 00448 
230 0. 00897 
3:40 0. 01796 
4.0 0.03598 
5:50 0. 04502 
6.0 070905 


ve 


* THE FOLLOWING BLOCK OF CODE INPUTS THE NODAL COORDINATES AND SPECIFIES 
* THE NUMBER OF ELEMENTS DESIRED BETWEEN NODES. 
SYSTEM 1 TYPE=CARTESIAN 


COORDINATES 
ENTRIES NODE K Z 

l 1. 49 4.48 
2 1.0 4.48 
3 0.0 4.48 
4 1.49 3.0 
5 1.0 3.0 
6 0. 0 2:70 
7 1.49 2. 24 
8 1.0 2728 
2 0.0 2.24 
10 1.49 0.0 
11 1.0 0.0 
12 0.0 0.0 

e 

LINE STRAIGHT 1 2 EL=3 

LINE STRAIGHT 2 3 EL=6 

LINE STRAIGHT 4 5 EL=3 

LINE STRAIGHT 5 6 EL=6 

LINE STRAIGHT 7 8 EL=3 

LINE STRAIGHT 8 9 EL=6 

LINE STRAIGHT 10 11 EL=3 

LINE STRAIGHT 11 12 EL=6 

Ye 

LINE STRAIGHT 1 4 EL=9 

LINE STRAIGHT 4 7 EL=5 

LINE STRAIGHT 7 10 EL=14 

LINE STRAIGHT 2 5 EL=9 

LINE STRAIGHT 5 8 EL=5 

LINE STRAIGHT 8 11 EL=14 
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LINE STRAIGHT 3 6 EL=9 
LINE STRAIGHT 6 9 EL=5 
LINE STRAIGHT 9 12 EL=14 
e 
* THE FOLLOWING CODE DEFINES MATERIAL PROPERTIES TO BE USED IN THE MODEL. 
THE FIRST MATERIAL SPECS CORRESPOND TO THE FIBER AND THE SECOND TO THE 
ALUMINUM MATRIX. EACH SET CONTAINS THE CORRESPONDING COMMANDS FOR 
* ELASTOPLASTIC ANALYSIS. THE "TREF" COMMAND MUST BE SET TO 530.0 FOR THE 
* QUENCH RUN AND 25.0 FOR THE SUBSEQUENT TENSILE DISPLACEMENT RUN. 
MATERIAL 1 THERMO-ELASTIC TREF=530.0 
0.0 450E9 0.2 4.3E-6 

1000.0 450E9 0.2 4.3E-6 
*MATERIAL 1 ELASTIC E=450E9 NU=0. 2 
* 


+ 


MATERIAL 2 THERMO-PLASTIC TREF=530. 0 
070 O 27616  656E6. 23. 0E-6 

1000. 0 emo 0.55. 27/6686 650656. 723, OE-6 
*MATERIAL 2 PLASTIC E=68. 3E9 NU=0. 33 Y=276E6 ET=656E6 
* THE FOLLOWING CODE SELECTS GROUPS OF ELEMENTS IN THE MODEL AND ASSIGNS 
fee) PREVIOUSLY DEFINED MATERIAL PROPERTIES 
EGROUP 1 TWODSOLID SUB=0 M=1 
O URE 5 6 9 8 EL1=6 EL2=5 NODES=4 
BE WE 6 9 12 11 EL1=6 EL2=14 NODES=4 
ye 


EGROUP 2 TWODSOLID SUB=0 M=2 


GSURF 1 2 5 4 EL1=3 EL2=9 NODES=4 
BE VRE 2 3 6 5 EL1=6 EL2=9 NODES=4 
GSURF 4 5 8 7 EL1=3 EL2=5 NODES=4 
GSURF 7 8 11 10 #£EL1=3 EL2=14 NODES=4 


se 


% THE FOLLOWING CODE CONTAINS THE BOUNDARY CONDITION SPECIFICATIONS. THE 
x FOLLOWING CONDITIONS ARE APPLICABLE FOR THE ALIGNED FIBER MODEL. 
* APPENDIX B CONTAINS A SECTION OF CODE ILLUSTRATING THE 
* APPROPRIATE MODIFICATIONS FOR THE STAGGERED FIBER MODEL. 
6011111 NODES 
1874 
BO 101111 NODES 
VOR 11 
34 TO 40 
207110111 NODES 
3 / 6 / 9 
SOTO 112 
* 
CONSTRAINT 
mo 2533 
E O LO 19 353 3 
2 1072 
4 2 10 2 
2 LO 2 
ime 0 2°70 64 2 10 2 
Je 
% THE FOLLOWING COMMANDS DEFINE THE TYPE OF LOADING TO BE APPLIED TO THE 
% MODEL. THESE COMMANDS ARE USED IN CONJUNCTION WITH THE TIMEFUNCTION 
* COMMANDS DISCUSSED PREVIOUSLY. IT IS IMPORTANT TO NOTE THAT IF THE 
APPLIED LOADING IS STATIC (I.E. TIME INDEPENDENT) THEN THE TIME VALUES 
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% IN THE TIME FUNCTION COMMAND REPRESENT LOAD STEPS. 
LOADS TEMPERATURE TREF=530 
LAO A TO 280 lene 
sv 
LOADS DISPLACEMENT 
36610 2 
sv 
* THE FOLLOWING COMMANDS ALLOW THE INPUT NODE AND ELEMENT CONFIGURATION 
* FOR THE MODEL TO BE DISPLAYED ON THE SCREEN AS LONG AS THE "WORKSTATION" 
* COMMAND IS COMMENTED OUT. 
FRAME HEADING=UPPER XSF=1.0 YSF=0.5 XFMAX=22.5 YFMAX=17.0 
MESH NODES=11 ELEMENT=1 
se 
ADINA 
END 
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APPENDIX B. ADINA MODIFICATIONS FOR STAGGERED FIBER 
MODEL 


The following modifications are used in the ADINA input file for the staggered fiber 
model analysis. Since many elements of the input file are unchanged from the aligned 


model, only required modifications are presented here. 


TIMEFUNCTION 4 

* DISPLACEMENTS FOR STAGGERED MODEL WITH RESIDUAL STRESSES 
0 0.0 

ea OL075 

500651 

. 00401 

. 00500 

700952 

103775 


Y Dd UY E O no 
ED DØD 
OOO 1 


* 


TIMEFUNCTION 5 

 DISPLACEMENTS FOR STAGGERED MODEL WITHOUT RESIDUAL STRESSES 
0.0 

0.00224 

. 00449 

. 00898 

301799 

202251 

, 04525 


DN Un SY NRO 
oO OC OO 
Soom © © 


* BOUNDARY CONDITIONS 
BO 101111 NODES 

7 
BO 110111 NODES 

eee? 9 / 12 

91 TO 115 
CONSTRAINT 
* TOP CELL SURFACE REMAINS PLANAR AND IS SLAVED TO NODE 12 

[meni O TO. 3 3 12 3 -1.0 

BE 10 TO 193 12 3 -1.0 
* BOTTOM CELL SURFACE REMAINS PLANAR 

NE MENI O To 11 3 12 3 1.0 

34 3 12 3 1.0 TO 40 3 12 3 1.0 
* AXIAL SIDEWALL CONSTRAINT 


NGS -1.0 
42 3 64 3 -1.0 
ĠIMJMOJT3 -1.0 
44 3 62 3 -1.0 
see 61. 3 -1.0 
5683 60-3 -1,0 
7.3 -1.0 
0058 3 -1.0 


dl 


4 
49 
50 
Sol 
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* 


U U L) U U 


5 
5 
5 
5 
5 


O $ U Qs y 


* LOAD COMMANDS 


LOADS DISPLACEMENT 


12437 = 17014 


sumi 
3 - 1 
9 pali 
JW 
pre 


O O O OO 
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APPENDIX C. DATA REDUCTION PROGRAMS 


The following programs were used to reduce the ADINA output file to data points 


to be plotted on stress-versus-strain curves. 


C THIS PROGRAM READS THE ADINA OUTPUT FILE AND COMBINES THE STRESS VALUES 
C WITH THE ELEMENTAL VOLUMES TO PRODUCE A VOLUME-AVERAGED VALUE OF AXIAL 
STRESS. 
C 

PROGRAM TEST 

CALL INIT 

CALL READ 

STOP 

END 
C 
C THIS SUBROUTINE ZEROS OUT ARRAYS AND READS THE VOLUME DATA FILE. 
C THE PARAMETER VARIABLES ARE: NUMBER OF ELEMENTS (NEL), NUMBER OF LOAD 
C STEPS (LC) AND NUMBER OF ELEMENTS IN THE FIBER (NFBR). 
C 


SUBROUTINE INIT 

PARAMETER (NEL=252,LC=6,NFBR=114) 

REAL SZ( 1: NEL,1:4),SZBAR( 1: LC) ,V( 1: NEL) ,FBR( 1: NFBR) 
COMMON // SZ,SZBAR,V,FBR 


DO 10 I=1,NEL 
V(1I)=0. 0 
DO 15 J=1,4 
SZ(1,J)90.0 
15 CONTINUE 
10 CONTINUE 


C 
DO 20 K=1,LC 
SZBAR(K)=0.0 
FBR(K)=0.0 

20 CONTINUE 

C 


OPEN (01,STATUSs'OLD') 

CORS b=! NEL 

READ (01,*,END=50) V(L) 
25 CONTINUE 


C 
50 PRINT *, 'EOF VOL. DAT' 
C 
RETURN 
END 
C 


C THIS SUBROUTINE READS THE FORMATTED ADINA OUTPUT FILE, AFTER THE HEADERS 
C HAVE BEEN REMOVED, AND CALCULATES THE AVERAGE Z STRESS IN EACH ELEMENT 
C BY AVERAGING THE FOUR STRESS VALUES AT THE GAUSS POINTS. THIS AVERAGE 
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MAMA 


25 
20 


50 


35 


STRESS IS MULTIPLIED BV ITS ELEMENTAL VOLUME AND SUMMED OVER THE ENTIRE 
VOLUME OF THE MODEL. THE OUTPUT IS TWO FILES: ONE CONTAINING THE 
COMPOSITE STRESS VALUES AT THE GIVEN LOAD LEVELS AND ONE CONTAINING 
THE STRESS LEVELS IN THE FIBER AT THE SAME LOAD LEVELS. 


SUBROUTINE READ 

PARAMETER (NEL=252,LC=6,VOL=4. 97 ,NFBR=114) 

REAL SZ( 1: NEL, 1: 4),SZBAR( 1: LC) ,V( 1: NEL) ,FBR(1: NFBR) 
CHARACTER*8 STATE 

COMMON // SZ,SZBAR,V,FBR 


OPEN (02, STATUSS OLD) 
OPEN (03,FILE='V30.DAT' ,STATUS='NEW') 
OPEN (04,FILE='V30_FBR. DAT' ,STATUSZ' NEW' ) 


DO 20 K=1,LC 

SUM=0. 0 

DO 25 L=1,NEL 

READ (02,*) IEL 

READ (02,21) INR,INS,STATE,SX,SY,SZ(L,1),SYZ,SVM,SMAX 
READ (02,*) SYS,SMIN 

READ (02,*) 

READ (02,21) INR,INS,STATE,SX,SY,SZ(L,2),SYZ,SVM,SMAX 
READ (02,7) SYS,SMIN 

READ (02,7) 

READ (02,21) INR,INS,STATE,SX,SY,SZ(L,3),SYZ,SVM,SMAX 
READ (02,*) SYS,SMIN 

READ (02,*) 

READ (02,21) INR,INS,STATE,SX,SY,SZ(L,4),SYZ,SVM,SMAX 
READ (02,7) SYS,SMIN 

READ (02,*,END=50) 

READ (02,%,END=50) 

FORMAT (14X,11,3X,11,5X,A8,5X,4(E12.5,2X),5X,2(E12.5,4X)) 


SBAR<((SZ(L,1)4SZ(L,2)4SZ(L,3)4SZ(L,4))/4. 0)%V(L) 
SUM=SUM+SBAR 


IF (L.EQ.NFBR) THEN 
FBR(K)=SUM/VOL 
ENDIF 

CONTINUE 
SZBAR(K)=SUM/VOL 
CONTINUE 


PRINT *, "EOF DETECTED” 


DO 35 N=1,LC 

WRITE (03,7) SZBAR(N) 
WRITE (04,%) FBR(N) 
CONTINUE 

RETURN 

END 
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C THIS PROGRAM CALCULATES ELEMENTAL VOLUMES BASED ON THE GEOMETRY OF THE 

Smet TE-ELEMENT MODEL. THE REQUIRED INPUTS ARE: 

NS; NUMBER OF ELEMENTS IN THE MODEL IN THE HORIZ. DIRECTION. 

NY: NUMBER OF ELEMENTS IN THE MODEL IN THE VERT. DIRECTION. 

NEA: NUMBER SG Rep DEER SB iEMENTS IN THE HORIZ. DIRECTION. 

NAN KOP EEK ELEMENTS TIN THE VERT. DIRECTION. 

NEM: NOMBEK OF ELEMENTS IN THE VERT. DIRECTION FROM THE 
MODEC'MIDEOINT TO THE FIBER END: 

NEL: NUMBER OF TOTAL ELEMENTS IN THE MODEL (NX x NY). 

AF: RADIUS OF FIBER DIVIDED BY NFX. 

XM: (RADIUS OF CELL - RADIUS OF FIBER) DIVIDED BY (NX-NFX). 

NE LENO SOG SE IUERK S07 LENGTH OF GELL) DIVIDED BY NFM. 

EZ NOS TESNO TVO CELE DIVIDED BY CNFY-NEM). 

Ke NOEL  LENGIH OF FIBER) DIVIDED BY I(NY=-NFY) 


C2) C2090) C2 RAÐ a a 


PROGRAM CVOL 
PARAMETER (NX=9,NY=27 ,NFX=6 ,NFY=18 , NFM=4 ,NEL=243, 
& XF=. 1666,XM=. 1633, YF1=. 19, YF2=. 16, YM=. 1644) 
REAL VOL( 1: NEL) ,V(1: NX, 1: NY) 
Ç 
C INITIALIZE THE MATRICES 
DO 10 I=1,NX 
DO 10 J=1,NY 
lowe V(1,J)20.0 
C 
DO 15 I=1,NEL 
15 VOL(1)=0. 0 
C 
C CONSTRUCT THE V MATRIX AND MAP IT INTO THE VOL MATRIX. 
C  *** TOP FIBER GROUP *** 
NID=1 
DO 20 L1=NFY,NFY-NFM+1,-1 
DO 21 K1=NFX,1,-1 
V(K1,L1)= ((XF*K1)%*2 ` (XF*(K1-1))%*2)*YF1/2. 
VOL(NID)=V(K1,L1) 
NID=NID+1 
21 CONTINUE 
20 CONTINUE 
C jek BOTTOM FIBER GROUP fik 
DO 200 L10=NFY-NFM,1,-1 
DO 210 K10=NFX,1,-1 
V(K10,L10)= ((XF*K10)**2 - (XF*(K10-1))**2)*YF2/2. 
VOL(NID)=V(K10,L10) 
NID=NID+1 
210 CONTINUE 
200 CONTINUE 
C ki TOP RIGHT MATRIX ELEMENT GROUP se 
DO 22 L2=NY,NFY+1,-1 
DO 23 K2=NX,NFX+1,-1 
V(K2,L2)=((XM*(K2-NEX)+XF*NFX)*%2- 
& ( XM*( K2-1-NFX)+XF*NFX)**2)*YM/2. 
VOL(NID)=V(K2,L2) 
NID=NID+1 
23 CONTINUE 
22 CONTINUE 
C  *** TOP LEFT MATRIX ELEMENT GROUP *** 
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DO 24 L3=NY,NFY+1,-1 
DO 25 K3=NFX,1,-1 
V(K3,L3)= ((XF*K3)**2 - (XF(K3-1))**2)*YM/2. 
VOL(NID)=V(K3,L3) 
NID=NID+1 
25 CONTINUE 
24 CONTINUE 
C der MIDDLE MATRIX ELEMENT GROUP **% 
DO 26 L4=NFY,NFY-NFM+1,-1 
DO 27 K4=NX,NFX+1,-1 
V(K4,L4)= ((XM*(K4-NFX)+XF*NEX)**2 - 
(XM*(K4-1-NFX)+XF*NFX)*%*2)*YF1/2. 
VOL(NID)=V(K4,L4) 
NID=NID+1 
27 CONTINUE 
26 CONTINUE 
C det BOTTOM MATRIX ELEMENT GROUP ton 
DO 260 L40=NFY-NFM,1,-1 
DO 270 K40=NX,NFX+1,-1 
V(K40,L40)= ((XM*(K40-NEX)+XF*NFX)**2 - 
5 (XM*(K40-1-NFX)+XE*NFX)**2)YF2/2. 
VOL(NID)=V(K40,L40) 
NID=NID+1 
270 CONTINUE 
260 CONTINUE 
C  *** CHECK TO ENSURE THE SUM OF THE ELEMENTAL VOLUMES MATCH THE siek 


C  *** CALCULATED MODEL VOLUME AND WRITE THE OUTPUT TO A FILE. AN 
VNUM= (C(C(XF*NEX)+(XM*(NX-NFX)))*%*2)*((YF 1*NFM)+ 
& CYF2*( NEY =-NEM) NNN 

C 
VCALC=0. 0 


DO 30 N=1,NEL 
VCALC=VCALC+VOL(N) 

30 CONTINUE 
PRINT *, 'VNUM="', VNUM 
PRINT ~, "VCALCSI SI CATE 
DONSSEMENMNEL 

35 WRITE (10,7) VOL(M) 


STOP 
END 
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